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Abstract 

In this article we compute and analyse the transition rates and duration of reactive trajectories of the 
stochastic 1-D Allen-Cahn equations for both the Freidlin-Wentzell regime (weak noise or temperature 
limit) and finite-amplitude white noise, as well as for small and large domain. We demonstrate that 
extremely rare reactive trajectories corresponding to direct transitions between two metastable states are 
efficiently computed using an algorithm called adaptive multilevel splitting. This algorithm is dedicated to 
the computation of rare events and is able to provide ensembles of reactive trajectories in a very efficient 
way. In the small noise limit, our numerical results are in agreement with large-deviation predictions such 
as instanton-like solutions, mean first passages and escape probabilities. We show that the duration of 
reactive trajectories follows a Gumbel distribution like for one degree of freedom systems. Moreover, the 
mean duration growths logarithmically with the inverse temperature. The prefactor given by the potential 
curvature grows exponentially with size. The main novelty of our work is that we also perform an analysis 
of reactive trajectories for large noises and large domains. In this case, we show that the position of the 
reactive front is essentially a random walk. This time, the mean duration grows linearly with the inverse 
temperature and quadratically with the size. Using a phenomenological description of the system, we are 
able to calculate the transition rate, although the dynamics is described by neither Freidlin-Wentzell or 
Eyring-Kramers type of results. Numerical results confirm our analysis. 
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1 Introduction 

This article is a numerical and theoretical analysis of reactive trajectories between two metastable states of 
the 1-D stochastic Allen-Cahn equation. Reactive trajectories are paths which directly connect two stable or 
metastable sets A and B (see Fig. [1]). Metastability is ubiquitous in physics and chemistry, for instance in solid 
state physics [T], physical chemistry [5] or molecular dynamics [3]. It has also been studied recently in new 
fields, for instance in fluid dynamics systems, like in 2-D Navier-Stokes equations j4], turbulent Von-Karman 
flows [5], turbulent convection [^, or wall turbulence [7]. In oceanic and atmospheric flows, well-known 
examples are the bimodal behavior of the Kuroshio paths [5], thermohaline circulation in the North-Atlantic 
liiin] or zonal-blocking transitions in the midlatitude atmosphere m ; in geophysics, one can also think of 
reversal of dynamos [T^]. Many of these systems can be modeled by stochastic partial differential equations 
(SPDEs). 
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In the last three decades, the stochastic Allen-Cahn equations and more generally reaction-diffusion 
equations have become central in the mathematical study of SPDEs [T3[Tl[T5l[T6l[l7l[T8]. They exhibit 
both relative simplicity and nontrivial behavior such as bimodality and effects of large dimensions. They are 
widely used for the study of phase-ordering kinetic (see e.g. mm]). In particular, the stochastic Allen-Cahn 
equations can be found in fluid mechanics as a convenient model of long-wavelength modulations in turbulent 
Couette flows (where it is termed Ginzburg-Landau equation) [7] , or in the context of quantum mechanics as 
a model of a double-well oscillator (see HD). In these systems, the metastable state comprise of a single phase 
in the whole domain. Meanwhile, during reactive trajectories the domain contains subdomains of different 
phases separated by fronts. 

For systems with one degree of freedom like stochastic differential equations (SDEs), solutions of the 
Fokker-Planck equations can often be obtained and analyzed analytically [Snilllllllj. As a consequence, all 
the statistics of reactive trajectories can be inferred. However, when the system of interest has too many 
degrees of freedom, an explicit use of such methods is often impossible. In such situations, several approaches 
can bring insight on the transitions. The discussion mostly depends on the noise amplitude. For relatively 
large noise, transitions between metastable states can be computed by direct numerical simulations. As the 
noise intensity decreases, these transitions occur with smaller transition rates and are therefore more difficult 
to observe using direct numerical simulations. In such a case, a possibility is to consider the small-noise limit 
theories such as the Freidlin-Wentzell [23] or the Eyring-Kramers [24l 123 EO Ell asymptotic theories for 
transition rates. The Freidlin-Wentzell large-deviation theory describes the most probable transition called 
instanton and the Eyring-Kramers expression gives the full formula, with the sub-exponential prefactor, for 
the mean first transition time or equivalently the transition rate. This formula holds strictly for gradient 
stochastic differential equations, but can be generalised (with modification) to a larger class of non-gradient 
system [28]. The case of SPDEs offers another level of difficulty: large-deviations results are very difficult 
to establish in general. For instance, it is not even known if a large-deviation principle exists for the Allen- 
Cahn equation in large dimension of space |18j . However, since the seminal work of Paris & Jona-Lasinio 
[13], the 1-D Allen-Cahn equation appears to be now well understood (see e.g. [IS] [E] [E] ), and follow the 
classical finite dimension phenomenology. Note that even for finite-dimensional systems, serious theoretical 
difficulties may also arise when the saddle points of the deterministic dynamics are degenerated |29j , or even 
worst, when the local attractors representing the metastable states are non-isolated [30]. Nevertheless, these 
large-deviation theories have been used with success on some SPDEs (see e.g. m) including the stochastic 
Allen-Cahn equations where finite-time instantons have also been numerically computed with great accuracy 
[32]. Another relevant approach of SPDEs is the direct computation of instantons [33l[34] using essentially a 
geometrical approach [35] . 

In the 1-D Allen-Cahn case, if the domain is not too small, instantons and transitions between metastable 
states take the form of fronts propagating from one side of the domain to the other. Previous results like 
the calculation of E et al. [32] as well as the small-noise limit theories in general do not take into account 
finite-amplitude noise effects. As will be discussed in this article, the range of validity of the Freidlin- 
Wentzell asymptotic regime is however extremely limited. The regime with noise too large to reach the 
Freidlin-Wentzell asymptotic, but still small enough for the event to be rare, occupies a very large area of 
the parameter space for partial differential equations. For symmetric potentials, a theoretical analysis of 
metastability for the I-D Allen-Cahn case can be tackled if the size of the domain is large compared to the 
front thickness [361I3Z1 EH ESlIig. In this case, the propagation of the fronts is a random walk. All things 
considered, the approaches are very diverse and leave whole areas of the parameter space unexplored. In 
particular, the limit of validity of Eyring-Kramers and Freidlin-Wentzell theories has never been considered 
precisely. We propose here to give a global description of the statistical behavior of I-D Allen-Cahn equation 
for arbitrary noise amplitude as well as finite-size domains. 

In order to do that, we consider a new numerical tool dedicated to the computation of rare events and 
in particular reactive trajectories. This algorithm is called adaptive multilevel splitting (AMS) and has been 
proposed rather recently by Cerou & Guyader m- It corresponds in fact to a modern adaptive version of 
the classical multilevel splitting algorithm introduced in the 50s by Kahn & Harris [42] and Rosenbluth & 
Rosenbluth [43] . Classical splitting algorithms have been very popular in the last few decades and have been 
applied in many areas of physics (chemistry, molecular dynamics and networks). This class of algorithms is 
now better understood from a mathematical point of view and have become an important area of research 
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Figure 1: Upper panel: a first passage trajectory (with mean duration called T in the text) is shown in 
black starting from a fixed point Xq and its final portion corresponding to the so-called reactive trajectory 
(with mean duration called r in the text) is shown in red. Lower panel: after generating JV i.i.d. trajectories 
starting from x, AMS is eliminating the trajectory with the worst excursion (here trajectory 1 corresponding 
to Qi) with respect to the reaction coordinate $ (see Appendix) and replaces it using information picked 
from a randomly chosen trajectories among the N — 1 “better” ones (here N = 3 and trajectory 2 is chosen 
by 50% chance). This is done by restarting a new trajectory from the initial condition (the red dot in the 
hgure). This process is repeated until all trajectories has reach the set B before A. 


in probability theory [44) . AMS appears to be very efficient and robust, it is moreover much less sensitive to 
tuning parameters than the classical versions and can in principle be applied to systems with a large number 
of degrees of freedom. In fact, the main constraint is to work with Markovian systems. The idea of the 
algorithm is rather simple to grasp: it is essentially a go-with-the-winner strategy (or successive iterations of 
the weakest link game) in the space of trajectories. At the start of the algorithm, one generates an ensemble 
of trajectories which are statistically independent and identically distributed. By successively eliminating the 
bad trajectories and generating new ones by splitting them from the good trajectories, one is able to obtain 
easily a whole set of trajectories connecting one metastable state to another. In principle, this can be done 
for any amplitude of the noise and possibly for more exotic colored noise as well. One is then able to perform 
statistics for rare events such as estimating the distribution of the duration of reactive trajectories. Although 
theoretical subtleties still remain regarding convergence properties [JS) US) 03 SH] ; AMS approach is likely 
to offer one of the most efficient tool for handling rare events computation. Until now this algorithm has 
only been applied to SDEs with a small number of degrees of freedom, and mostly for test and convergence 
analysis. One of the main aim of this study is demonstrate the practicality and usefulness of this algorithm 
in high dimensional and spatially extended systems as well. Moreover, as will be shown, the numerical results 
will provide a much broader view on the theory of reactive trajectories. A sketch of the algorithm is shown 
in FiglU A detailed discussion is also given in the Appendices IA.2I and IA.3I together with a pseudo-code 
description. 

This article shows several results. We are able to obtain a complete theoretical and numerical statistical 
description of the reactive trajectories for the 1-D stochastic Allen-Cahn equation in the parameter plane 
spanned by the inverse temperature P and the length L of the domain. Those are the only control parameters 
of the Allen-Cahn equation once properly rescaled. We also discuss the type of reactive trajectories, the 
transition rate, the mean first passage time and the duration of reactive trajectories as a function of /3 and L. 
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We show that only a small portion of the parameter plane for large j3 and small L leads to ensemble of reactive 
trajectories that concentrate close to the most probable one, as could be expected from Freidlin-Wentzell 
theory and instanton phenomenology. We investigate the transition rate and the mean first-passage time and 
recover the classical Eyring-Kramer law in this narrow band of the parameter plane. We show in particular 
that the prefactor has an exponential dependence in the size. We are also able to distinguish the limit of 
validity of this law. In some parts of the parameter space, for finite-amplitude noise, reactive trajectories are 
close to instantons from a geometric point of vue, but their temporal evolution is not close to the instanton 
one, by contrast to the Freidlin-Wentzell regime case. In particular, they are characterized by a random-walk 
propagation of their fronts. For a very large size, the system presents the properties of a potential with a 
extremely flat saddle. This may be related to the case of non-quadratic saddles for which classical results 
breakdown pUl H5] . Note that this can also be studied in a case where the global minima of the potential do 
not have the same height [13] • The dynamics is described by an effective potential whose curvature at the 
saddle in exponentially small in L. We can write and solve a coupled master/Fokker-Plank equations that 
describes the position of the front. From this, we can calculate the transition rate and show that its prefactor 
is not exponential in L, as would be expected from Eyring-Kramers expression with an exponentially small 
curvature, but rather polynomial. However, the transition rate retains its large deviations properties in /3. 
This analytical result is verified by numerical computations using AMS. We are also able to investigate the 
statistics of the reactive trajectory durations as well. We find two distinct regimes of large curvature and 
small curvature at the saddle where the duration scales like exp{L/^/2) ln/3 and L^j5 respectively. Moreover, 
in the large curvature regime, the distribution of durations follows a Gumbel distribution like for one degree 
of freedom systems [50] . As far as we know, these are new results. 

The structure of this work is as follow. We first introduce the I-D stochastic Allen-Cahn and its de¬ 
terministic solutions in Section We then perform in Section [3] a theoretical analysis of metastability in 
Allen-Cahn equation using small noise theories such as Freidlin-Wentzell large deviation theory and Eyring- 
Kramers law. Section |3| extensively makes use of the AMS algorithm for computing numerically ensemble of 
reactive trajectories in various regimes of the Allen-Cahn equation. We will be particularly interested in the 
comparison between theory and numerics in the range of parameters where both are valid. We discuss these 
results in the conclusion {%U- 

2 The stochastic Allen-Cahn equation 

We present here the I-D stochastic Allen-Cahn equation with a definition of the control parameters and some 
properties of the stationary solutions of the associated deterministic Allen-Cahn equation. 

2.1 Introduction 

The nondimensional I-D stochastic Allen-Cahn equation with Dirichlet boundary conditions reads 

dtA = d’l.A + [A — A^) -I- ^(0) = A{L) = 0, (I) 

where 77 is a white noise in space and time, < rj{x,t)ri{x',t') >= 6{t — t')6{x —x') |5T1|52]. It has several names: 
it is often called the Allen-Cahn equation in mathematics and solid state physics [1] and the Chaffee-Infante 
in climate science. In non-linear physics, it is one of the several instances of Cinzburg-Landau equation [53]. 
In the study of the deterministic front propagation in biology, it is a particular case of the Fisher-Kolmogorov 
equation [54] . After change of temporal, spatial and amplitude scales, there are only two free parameters, 
the nondimensional inverse temperature /3 and the nondimensional size of the domain L. With this scaling 
choice, the size of the boundary layers, fronts and coherent length is of order one. We recall an important 
property of the Allen-Cahn equation which is essential for theoretical considerations: it is a gradient system 
(for the scalar product) with respect to the potential V 

^ ^ ^ i ^ T + 
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This formulation allows one to have an easier understanding of the phase space structure and dynamics of 
the Allen-Cahn equation. 

2.2 Stationary solutions 

Stationary solutions of the deterministic Allen-Cahn equation solve the Sturm-Liouville problem (see Paris 
& Jona-Lasinio [S] §7 and references therein) 

dlA=-^, A(0) = A(L) = 0, (3) 

where U = A^/2 — A^/A. Depending on L, several solutions may exist, the stable ones play the role of the 
metastable states in the stochastic Allen-Cahn equations whereas saddle points indicate possible pathways 
between them. Equation ([3]) can be solved using the analogy of a pendulum of period L in a potential U 
where x plays the role of time. Provided L is not too small, there are two stable solutions denoted and 
one unstable solution Aq (Fig. [2]) 

Ao = 0,Aq =±1 + “boundary-layer corrections”. (4) 

The two stable solutions correspond to the global minima of V (Eq. ([^l ) ■ Due to the existence of Dirichlet 
boundary conditions at a: = 0, L, boundary-layer for should also be considered. As L is increased, pairs of 
unstable solutions denoted (—A„, A„) successively appear. It has been demonstrated in [13] that the critical 
values of L for which this happens are 

L = {n + l)7r (5) 

These solutions have exactly n fronts (called kinks in field theory [T]) and n zeros in the domain ]0, L\. Using 
the dynamical analogy of the pendulum, it means that the solution are periodic with period 2L/{n + 1). Note 
that the Allen-Cahn equation without boundary conditions and the diffusive part d^A has a whole family 
of stationary solutions solving A^ = 1. These solutions can arbitrarily jump from the values -1-1 and -1 at 
any points of the domain. The diffusive term together with the boundary conditions select some of these 
solutions. 

Figure [3] shows three of the solutions of (j3|) for L = 30 where 1 -|- 2 x 8 solutions coexists due to the 
symmetry A —>■ — A (see (O). The figure shows the stable solution Aj and two saddle solutions Ai (one 
front) and A2 (two fronts). The size of the fronts and boundary layers is 0(1) as expected and the solutions 
take values close to ±1 away from the boundary layers. 

It is possible to give an analytical expression of the saddle Ai locally near x = L/2 for large L. We 
denote A^ the solution which corresponds to a solution of 9^A = —dU/dA on K. with limit conditions 
lima;_,.±oo Ai(x) = ±1. It takes the exact form 

Ar{x)=tanh(±^y (6) 

This is indeed the classical kink formula HI [36]. The term ± comes from the symmetry of the problem 
A —>■ —A, it is -|-1 for an ascending front or kink and —1 for a descending front or anti-kink. A good 
numerical match is found in Fig. [5] between this expression and the exact solution Ai away from the external 
boundary layers. Note that for more complex potentials, Taylor expansion of the potential would yield local 
expressions as well (see [13] § B). 

The analysis presented by Faris & Jona-Lasinio m indicates how many fixed point the deterministic 
system has and whether they are stable. In the case of Dirichlet boundary conditions, it does not give the 
eigenvalues of the Hessian of V at said fixed points: these eigenvalues indicates how stable or unstable the 
hxed points are. Using analytical description of the type of Eq. ([5]) we are able to give an approximation of 
the one negative eigenvalue of the Hessian of V at Ai (the most important unstable fixed point for reactive 
trajectories). We describe the result in section 13.2.2 1 and detail the derivation in appendix lA. 1 1 In particular, 
we are able to show that this eigenvalue converges toward 0 exponentially with a rate exp(—L/-\/2). This 
indicates that the potential V becomes extremely flat at the saddle Ai. This flatness is the main cause of the 
diffusive regime described in section l3Al 
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Figure 2: Example of stationary solutions of the 1-D deterministic Allen-Cahn equation in a domain of size 
L = 30. The solution Ai with one front is compared to the analytical solution of the front Eq. (jO)). Equation 
m is solved by finite-difference using a Newton method with initial condition sin7r(n -|- l)a;. 


3 Theoretical methods for metastable systems 

We present here the theoretical approaches we use to analyse metastability. We also complement these by 
specific discussions and results whenever it is necessary. First, we discuss the Freidlin-Wentzell large-deviation 
approach in subsection l3.ll then Eyring-Kramers formula and the distribution of reactive trajectory durations 
in subsection l3.2l then the front dynamics description in subsection l3.3l and summarize these results in a phase 
diagram in subsection 13.41 

3.1 Large deviations, Freidlin-Wentzell theory, and Instantons 

We consider a gradient dynamics in a potential V with a weak noise of amplitude •^2//?. We consider the 
transition from the basin of attraction of one local minima uq of V to another one ui. We assume that the 
two basins of attraction touch each other at at least one saddle point. Generically the transition rate fj, from 
the basin of attraction uq to the one of ui follows a large deviation result 

lim - - y(uo), (7) 

/3->-oo p 

where V{us) is the minimum of V at any saddle points connecting the two basin of attraction of uq and ui. 
When the lowest saddle is unique, the reactive trajectories concentrate close to a single path called instanton- 
anti-instanton. This phenomenology can be proven within the Freidlin-Wentzell theory [53]. As far as the 
ID-Allen-Cahn equation is concerned, the proof of large deviation principles has first been established by 
Faris & Jona-Lasinio m- In this section, we explain heuristically those results, and describe the classical 
results about the saddle points and instantons of the ID-Allen-Cahn equation. 

We consider a partial differential equation 


du = F{u)dt + 



( 8 ) 


where u{t) = u(x, t), x G u{t = 0) = uo(x) is an initial condition and Wt is a (space-time) Wiener process. 
We call V = V(uo,mi,t) the set of trajectories which start from uq and end at ui at time r. A path integral 
approach (Initiated by Onsager and Machlup [55]) allows to formally express the transition probability from 
the state uq at time t = 0 to the state ui after time t = t as 


P{ui,t;uo,0) 




V[u], 


(9) 
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where the action is 


2 




du 

dt 


-Fiu) 


dt. 


In the low noise limit, the transition rate /r between attractors Uq and Ui is defined as P{ui,t; Mq, 0) ^ 

^T. From the path integral ([3]), taking first the limit /3 to infinity through a saddle-point approximation, and 
then the limit r —>■ oo yields the large-deviation result 


lim — — InP = - inf £[m, oo], V = {u(t)\ lim u{t) = uq , u{0) = Us} 
/ 3 —>00 p 2 ugv t—> —00 


( 10 ) 


where Us is a saddle point at the boundary of the uq basin of attraction. Remarkably, Freidlin and Wentzell 
have shown that this saddle approximation is indeed valid provided uq is a fixed isolated attracting point of 
the deterministic system. Moreover, for gradient systems where F{u) = — the action minimization is easy. 
For any Ui in the basin of attraction of uq, we call a relaxation path from Ui to uq a solution of 


du 5V 

dt Su 


( 11 ) 


with initial condition u(t = 0) = Ui. We call a fluctuation path a time reversed relaxation path |56] . Then the 
unique minimizer of the infinite time action from Ui to uq is a relaxation path, the action minima is then equal 
to zero ; and the unique minimizer of the infinite time action from uq to Ui is a fluctuation path, and the action 
minima is then equal to 2V {ui) — 2V (uq) [53]. We call an instanton a solution of (fTTl) with limt_,,_oo u{t) = Uq 
and limt_,,oo u{t) = Ug (the instanton is then defined up to an arbitrary time shift). An anti-instanton is a 
time reversed instanton. Then it is easy to understand that inf^igv 00 ] = 2V(us) — 2V{uq), and that 
formula m leads to the result 0. Because the saddle point Ug is at the boundary of the basin of attraction 
of Uq, we stress however that there is not strictly speaking a relaxation path from t = 0 at Ug to mq- As 
a consequence, no minimiser exist for Eq. (uni). Any path that approximates an instanton from uq to Ug 
followed by an anti-instanton from Ug to iti has an action that is close to 2V(ug) — 2V(uo). We discuss this 
further with the Allen-Cahn example below. 


We now apply these results to the Allen-Cahn equation using ug = and ui = Ag and for different 
saddles. For instance, for Ug = Ag = 0, Fig. |3| is a space-time representation of the succession of an 
approximate instanton-anti-instanton trajectory. One obtains a global flip of the solution. We refer such a 
solution as a “flip”. We note that the duration of the trajectory close to the saddle point (the green part on 
Fig. [3|) is chosen arbitrarily and that the minimum of the infinite time action is obtained in the limit where 
this duration is infinite. When the saddles are A„ for n > 1, the instantons-anti-instantons correspond to 
some fronts propagating either from the boundaries (Fig. [SJd) or in the bulk of the domain (Fig. [3]:). We will 
term “front” trajectories such solutions. The sharp transition which can be seen in Fig. [3]d,c corresponds in 
fact to a timescale much smaller than the motion near the saddle and not to some numerical artifacts. For 
these cases too, the time spent by this approximate instanton-anti-instanton trajectory close to the saddle 
point is arbitrary. From the formal interpretation of the large deviation through path integrals, as described 
above, we expect the distribution of reactive trajectory to concentrate close to approximate instanton-anti- 
instanton trajectories. For the actual dynamics with small but finite noise, the time spent close to the saddle 
point is however not arbitrary. This is a random variable, with an average value that scales like ln(,5) in the 
limit of small noise 1//3, as discussed further in section [3.3.11 

From the previous discussion, as several saddles exist, several instanton-anti-instanton trajectories can 
join Ug = Aq and ui = Aq . In the low noise limit, the relevant one is the one that minimizes V{ug) — V{ug). 
For large L, the potential difference AI 4 , = E(A„) — V{A^) for n > 1 can be investigated analytically. The 
idea is simply to notice that A„ and A^ mostly differ in the neighborhood of the fronts provided L is large 
enough so that these fronts are well separated. In practice, this requires that L ^ n\/2. Let us look at 
the case n = 1, the solution (see eqU]) is a good approximation of the saddle so that in this case using 
|AJ| ~ 1, one can write 

AEi + \{d,Ar f + dx. 
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Figure 3: Space-time representations of relaxation-fluctuations paths for three different saddles: a): Aq = 0 
(L = 6), b): Ai {L = 10), c): 4 I 2 {L — 30). Only a) and b) can be instantons. 



Figure 4: Numerical computation of the potential difference as a function of L for Ag = 0 et = 0. 

After simple algebra, the above integral gives AFi ~ The case of several fronts follows the same line of 
reasoning, each front gives a similar contribution in the potential difference provided L is large and Anally 

AK —n^^,n>l. (12) 

The conclusion is that minimizers can only be instanton-anti-instantons going either through Aq or Ai where 
the potential difference is the smallest. A more precise analysis shows that for L < 2tt, the saddle Ai does 
not exist. It appears through a pitchfork bifurcation for L = 2tt. For L > 2tt, AVi < AVq as illustrated 
on figure m We thus conclude that for L < 2tt only the flip instanton exists, while for L > 2tt the relevant 
instanton-anti-instanton is the one front transition. 

3.2 Transitions rates, mean first passage times and Eyring—Kramer s’law 

3.2.1 General discussion 

We consider the stochastic differential equation (or partial differential equation) ([8]) , assuming in this section 
that the force is gradient: F{u) = —W{u). We want to compute, asymptotically for small noise (/3 —>■ 00 ) the 
mean transition time to go from one potential minimum say Ug to another ui assuming that the two basins of 
attraction are connected through a saddle point Ug- If we have two local minima only, any trajectory starting 
from Ug will eventually hit a small neighborhood of ui, after having possibly wandered around ug for a very 
long time. We recall that the time for a trajectory to reach ui, is a random variable T. The mean first 
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passage time T = Eq [7^ is equal to the inverse of the transition rate T = l//r, asymptotically for small noise. 
The trajectory between t = 0 and t = T are called a first passage trajectory. 

Provided the saddle Us is non degenerate, the Hessian at Us has a unique negative eigenvalue \s(us) 
(corresponding to the dim-1 unstable manifold). It is a classical result that asymptotically for small noise 
(see [241 EH ESI Ez]), the mean first passage time is equal to 


T = 


27r 

|As(Us)| 


|detV2p(u^)| 

detV2p(uo) 


exp {(3{V{us)-V (uo))). 


(13) 


Geometrically, the negative eigenvalue As of the Hessian corresponds to the lowest curvature of the potential 
V at the saddle and controls the time scale of the dynamics around it. We note that this Eyring-Kramers 
formula generalize, in this gradient case, the large deviation result discussed in the previous section. The 
distribution of T is in general exponential, i.e. is Pr(T ~ s) = ^exp(—s/T) (see (IT^ l [5^ . 

A difficulty is to give a meaning of these results in the context of SPDEs when u now belongs to a 
functional space. For instance, it is easily seen for the ID Allen-Cahn equation that the product of the 
Hessian eigenvalues diverges. However, it can be shown that in the case of 1-D Allen-Cahn equation, the 
ratio which appears in eq. (USD can be well-defined (see [HI [161 [17] ). Another level of difficulty can appear 
however in dimension larger than 2 at the level of the large deviations. First, it appears that Allen-Cahn 
equations are in fact ill-defined when the spatial dimension is strictly larger than one (which is not the case, 
we have considered in this work) |141l57j . One has to renormalize the equation properly so that large-deviation 
results can then be obtained in dimension 2 and 3 which are very difficult to demonstrate |18] . In fact, the 
question is largely open in dimension larger than 3. At the level of the prefactor of Eyring-Kramer law nothing 
is known even in dimension 2. Only the one-dimensional case is well-understood [ii[iii[n]. 


3.2.2 Application to the Allen Cahn equation 

We now discuss the behaviour of the prefactor 


2^ l \nZ2^Ms)\ . 

vixOTVin“iA.(^o)i ^ ^ 

of the mean first passage time. A general property of such systems is that As(A 5 ) goes to zero in the infinite 
size limit while the rest of the prefactor is unchanged. We therefore concentrate our effort on that eigenvalue 
As (As) < 0: it controls the size dependence of the mean first passage time. 

This eigenvalue As is the result of the eigenvalue problem 


= -^2$, + (SA^ o - 1)$, = 


(15) 


arising from the Hessian of V, with boundary conditions 4>i(0) = ^i(L) = 0. Since this operator is self 
adjoint, a good approximation is given using the min-max Rayleigh-Ritz theorem: the Rayleigh quotient 
gives an upper bound on As. This upper bound can be as close as possible to As. A analytical description 
can thus give a very good approximation of As and in particular of its decay rate. We detail the calculation 
in appendix lA.II At first order in exp(—L/\/2), the approximation reads 

/~-24e"^. (16) 


The logarithm of the absolute value of As (calculated numerically) and I are displayed in figure [5] (a). The 
ratio / gives a very precise upper bound, the ratio between As and / is approximately 2, independently of 
L. This shows that the factor in the exponential is exactly the decay rate of As. Meanwhile, 24 gives a 
reasonable order of magnitude of the factor in front of the exponential. Note that a scaling in exp(—L) had 
been calculated for the rate of evolution along the slow manifold of infinite size Allen-Cahn systems when 
two fronts are separated by a distance L [5^ . 

Numerical computation of the whole spectrum of the Hessian of V also allows us to compute the rest of 
the prefactor (Fig. [S] (b)). We can see that it quickly converges toward an asymptotic value independent of 
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Figure 5: (a) : logarithm of the eigenvalue As and its approximation / as a function of the size L. (b) : 
logarithm of the rest of the prefactor as a function of size. 


L. This confirms that the dependence on L of the prefactor Eq. (fHl) comes solely from Ag. As a consequence, 
in the range of parameters in which the Eyring-Kramers law is valid, the mean first passage time will grow 
like exp(L/(2-\/2)) with the domain size L. 

3.2.3 Duration of trajectories in one degree of freedom systems 

The duration s of a reactive trajectory is defined as the time elapsed between the instant where A leaves C and 
the instant where A touches dB (see Fig. [T]). For systems with only one degree of freedom, the distribution 
of durations of reactive trajectories is known in the limit /3 —>■ oo (see |50jl: it is a Gumbel distribution with 
parameters a and b 


1 2 2 

/ N 1 / / 7 2 ^ 

Pa, bis) = -exp [-Sab - exp{-Sab)) ,T = mean = & + 7 a,var = ——, 
a D 


(17) 


with Sab = and 7 ~ 0.6 is the Euler constant. The normalized distribution of the duration pis) = pi^^) 
reads 




(18) 


The average r can be calculated as a function of /3, |As| (the opposite of the second derivative of the potential 
at the saddle). In fact, it has been demonstrated by Cerou et al. that in the case of a one dimensional 
gradient system with a non degenerate saddle, we have in the /3 —>■ 00 limit: r = (ln(/3) + ln(|As|) + C'o)/|As|. 
In this result, the constant Cq depends on the properties of the starting and arrival points |50] . In the one 
degree of freedom case, one can also show that the variance is independent of /3 and equal to 7 r/(-\/ 6 |As|). 
Performing such an analysis in a many degrees of freedom system such as the Allen-Cahn equation is beyond 
the scope of this article, mainly because it requires a precise knowledge of the committor function (see § IA.2I 
as well as [iniiiziisi] and references within) in the neighbourhood of the saddle. This function of A gives the 
probability that a realisation of the dynamics starting from a given A reaches dB before dA. 

A transposition of the one degree of freedom result on the movement of the front on the saddle would lead 
us to expect that r is affine in ln(/3) 

ln(/3) + ln(|As|) + Ci 


I As 


(19) 


with I As I oc e as approximated by I Eq. (fTHll . 


3.3 Transition rates and reactive trajectories beyond Freidlin—Wentzell and Eyring- 
Kramers regimes 

For large domain size L, fronts form in the bulk of the domain and have an effective interaction with the 
boundary and between themselves that decays exponentially with L. For small but finite noise, in the limit 
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of large domains, we expect the front dynamics to be mainly diffusive. In the infinite domain limit, it has 
actually been proven that the motion of the front is brownian with a diffusion coefficient D = 8/(3/3) [37] , 

We expect this diffusive front dynamics to be also relevant for reactive trajectories. This will be illustrated 
also by numerical computation in sectionlH This phenomenology contrasts with the classical Freidlin-Wentzell 
picture for which one expects the transition probability to be dominated by an action minimum contribu¬ 
tion and small fluctuations around it. Actually, as we explain in the following section, the saddle point 
approximation leading to a large-deviation result is valid only in the parameter range /3 ^ exp(L/-\/2)/T^. 
Similarly, the Eyring-Kramers formula for transition rates and mean first exit time is valid only in the same 
parameter range. Those remarks stress that the range of validity of the Freidlin-Wentzell regime and of the 
Eyring-Kramers formulas are quite limited for large domain sizes. 

In this section we show that in the parameter range exp(L/'\/2)/T^ ^ /3 ^ I/AE, even if the saddle 
point approximation fails, a large deviation result and an asymptotic formula for the transition rate can still 
be obtained through a specific approach. It is not a consequence of the Freidlin-Wentzell Principle of large 
deviations and the reactive trajectories are notably different from instantons. The transition rate differs from 
a Kramers-Eyring formula. We also discuss the typical duration of reactive trajectories and their distribution, 
which are markedly different from the classical ones in the Freidlin-Wentzell regime. 


3.3.1 Limits of Preidlin—Wentzell and Eyring-Kramers regimes 


In this section we discuss the range of parameters for which the Freidlin-Wentzell and Eyring-Kramers regimes 
are expected to be valid, first for a one dimensional stochastic differential equation, and then applied to the 
Allen-Cahn equation. 

We consider the one dimensional dynamics dC,/dt = —dVe/dC, -I- For this case, it is very classical 

to write quantities of interest in terms of explicit integrals. Those integrals can be evaluated asymptotically 
for small noise using a saddle point approximation |51) . For instance the Eyring-Kramers formula for the first 
passage time can be readily obtained this way. This derivation, or any other derivation, stress the importance 
of the dynamics close to the saddle point. The dynamics close to the saddle point can be described by 
dQ/dt = — AsC + V2Dr]{t). We consider this equation with initial condition C(0) = 0. What is the average 
of the first hitting time r for which ICKt) = L/27 What is a typical behavior of the trajectories leading 
to ICI(t) = A/2 ; is it of a diffusive nature (the term —As is negligible), or is it rather dominated by the 
instability (—As dominates, except for very small 0? A simple dimensional analysis is enough to prove that 
the average of the hrst hitting time is E(t) = f{D/\Xs\L'^)/\Xs\, where / is a non dimensional function. In 
the small diffusion limit D <C |As|L^, one recovers the behavior corresponding to the validity regime of the 
Eyring-Kramers relation. Then the distribution of r is a Gumbel distribution, as discussed in section 13.2.31 
with E(r) = C 2 ln(I3/|As|L^)/|As|, where C 2 is a non dimensional constant. In this case, after an initial 
short diffusive behavior, the trajectories follow an exponential escape. This regime corresponds also to the 
large deviation regime, where trajectories concentrate close to instanton-anti-instanton trajectories. In the 
opposite limit D ^ |As|L^, the dynamics is dominated by diffusion. Then E(t) = CsL^/D, with C 3 another 
nondimensional constant. In this regime the trajectory are diffusive and they do not have any concentration 
property. 

In order to apply this discussion to the Allen-Cahn equation, we first note that in the case of large 
domains, the dynamics can be characterized by the motion of fronts. A single front motion is described by 
an effective diffusive dynamics in an effective potential d(^/dt = —dVe/dC, + ^/2Dr]{t), with D = 8/3/3. As 
discussed in section 13.2.21 the effective potential I 4 decays exponentially fast away from the boundary, then 
As = Xo{L) exp{—L/y/2). The factor Aq contains the exact value of the slowly varying factor of this eigenvalue 
(see ? 13.2.21) . From the previous discussion, we conclude that the large-deviation and Eyring-Kramers regime 
are observed for 


/3 » f3*{L) 


A^IAol 


( 20 ) 


where is the length scale of the front, see for instance Eq. (HI). This is merely the explicit writing of the 
condition D <C |As|L^. We can alternatively write this condition L <C L*{j3), with L* = /3*“^(/3) the unique 
root of this equation. 
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3.3.2 Mean first passage times and transition rates for large domains durations of reactive 
trajectories 

We now consider large domains L ^ 1 and the range of parameters exp{L j^/2) jL? /3 ^ l/Ai^ for which 
the Eyring-Kramers formula and the Freidlin-Wentzell phenomenology is not valid anymore. In that range 
of parameter, the dynamics is dominated by the nucleation of one or several fronts and their diffusive motion. 
In this subsection we discuss the simplest case, when a transition between the two attractors is due to the 
nucleation of a single front at one boundary, its diffusion throughout the domain until it reaches the other 
boundary to form the second attractor. The generalization of this discussion to the dynamics of multiple 
fronts would be rather easy. In section 13.3.41 we will discuss the range of parameters for which the one 
front transitions are the typical ones. In this section, we propose a phenomenological theory of one front 
transitions, compute the transition rate (or mean first exit time) and discuss the statistics of durations of 
reactive trajectories. 

We call A and B the two attractors, and resp. Pg the probability to observe them. We make several 
assumptions which are very natural. First we assume that when the system is in the attractor A (resp. B) 
the nucleation of a front close to either one or the other boundary has a Poisson statistics characterized by a 
nucleation rate exp(—/3Ay)/x, with y a time scale that has a finite limit in the limit of small fUAV. This is 
natural as the nucleation is an activation process. Moreover it is natural to assume that x has a finite limit 
for L going to infinity. Secondly, we assume that in the asymptotic large L limit, the dynamics of the front is 
the same as in an infinite domain, studied for instance by Brassesco et al. 137]. The front then diffuses like 
a Brownian particle with diffusion coefficient D = 8/(3/3). Let us call P{x,t) the density of the probability 
to observe a front at distance x (resp. x — L) from its nucleation boundary if it has been nucleated from the 
attractor A (resp. B). We note that a front can be nucleated either at the right or the left boundary, from 
either the attractor A or the attractor B. However, because of the symmetry of the problem, and because we 
assume that their is a single front at a time, the probability P{x, t) defined that way is enough to characterize 
the state of the system. 

A key point is to relate the boundary condition for P{x,t) with P _4 and Pg. The process at the boundary 
involves the creation and the destruction of a front. In order to model these phenomena, we denote v the 
rate per unit of length for a front in the vicinity of the position a; = 0 disappears and form the attractor A. 
By symmetry, this velocity v is also the rate of disappearance per unit of length of a front close to a: = L to 
form attractor B. It is natural to assume that v has a finite limit when L goes to infinity, as the effect of the 
second boundary becomes negligible. The phenomenological equations then read 


dPA 

dt 

dPs 

dt 


exp(-/3Al/) 

X 

exp(—/3A1/) 
X 


Pa + vP[t)A ), 
Pb + vP{L, t) , 


coupled to the Fokker-Plank equation 

^ _ g^P 

dt dx"^ ’ 

with boundary conditions at a: = 0 and x = L 


D^(0.t) = + „p,o,,) .„d 

dx X 

dP, , exp(—,5Ay) 

= - -Pb - vP{L,t) ■ 

dx X 


( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 


The equilibrium distribution should be the one corresponding to the Gibbs distribution for the Allen- 
Cahn model. In the limit of infinite domain sizes, we expect, at equilibrium Pa = Pb ^ 1/2, and P{x) = 
exp{—f]AV)/li, where li is a length of order of the front size that could be computed easily from the Gibbs 
distribution using a saddle point approximation. Looking at the equilibrium of the phenomenological model 
()21II25|1 . we conclude that the relation li = 2vx must hold. We note that our phenomenological model verifies 
detailed balance with respect to its invariant measure. 
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While the value of D has been computed from the Allen-Cahn equation, the values of v and x remains 
phenomenological parameters so far. It would be a very interesting problem to compute them from the Allen- 
Cahn equation (as explained above they are related through li = 2vx)- 


We first compute the flux j = —DdP/dx for a steady solution of this phenomenological model for which we 
impose Pa = 1 and Pg = 0 (we consider equations (I23I24I25I) with P 4 = 1 and Pg = 0). From the steadiness 
of the process and Eq. (ESI), we have that P{x) = —jx/D + P(0). The flux j and P(0) are determined from 
the boundary conditions Eq. (IMl) and Eq. (1251) . We obtain 


J = 


D 


ui{vL + 2D) 


exp(—/3AI/) 


16 

L>1 piiL 


exp(-/3AE), 


where we have used for the left hand side D = 8/3/3 and 2vx = h- 

Such a steady solution is established after a diffusion characteristic time td = L'^/D = 3PL‘^/8. Whenever 
this diffusion time is much smaller than the nucleation time Xn = X exp(/3AE), it is easily understood that the 
evolution of P _4 and Pg is described by an effective slow equation. This slow equation can be easily derived. 
Hence for 1 <C T \/x!P exp(/3AE/2), we obtain 

^ = -M(i/ 4 -^^B)and (26) 

at 

^ - Pa) , (27) 

where ^ takes the asymptotic value of the current j computed above 

M = exp(-/3AE). (28) 


As a conclusion, whenever 1 ^ L <C \/xlP exp(,5Ay/2), the dynamics of the transition from the attractor A 
to the attractor B is an effective Poisson process, with transition rate //. The mean first exit time from the 
basin of attraction of one of the attractors is then l//i. 

We do not discuss the parameter range L y^x//3exp(/3AE), because the hypothesis of a one front 
solution is no more valid in that case. This limit is studied in section [3.3.41 We discuss in next section the 
range of parameters for which the one front states are actually the typical ones for reactive trajectories. 

We have discussed in this section a phenomenological model of the kinetics of a one front reactive trajectory. 
However, it should be noted that this model can be easily extended to a kinetic model including the possibility 
to have n fronts at the same time. 


3.3.3 Duration of reactive trajectories for large domains 

We now discuss the duration of reactive trajectories and come back to the range of validity of our hypothesis 
that only one front forms, always for large domain sizes. 

As it is clear from the phenomenological model described in the previous section, the kinetics of the 
transitions are governed by diffusion in the large size limit. Eor random walk in dimension 1, it can be 
demonstrated that the average duration of reactive trajectories is proportional to the square of the size of the 
domain and inversely proportional to the diffusion coefficient [50]. We then have 

r = ( 74 / 3 X 2 , (29) 

where c is a non-dimensional constant. Note that in the case of a free random walk with one degree of 
freedom, the expected pdf of trajectory durations is not a Gumbel distribution any more. Moreover, the 
variance depends on /3. In fact, it is proportional to the average duration var = t-\/2/5 [50] . 

The one front case will remain the typical transition state of the system only if the average duration of 
reactive trajectories r = ( 74 / 3 P^ is much smaller than the nucleation time Xn = X6xp(/3AH). We conclude 
that this is true whenever 

V^exp(^^) . 
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This criterion turns out to be the same as the one derived for the validity of effective Poisson process for 
describing transitions from one state to another. The same criteria will be also obtained in the next section 
using a purely equilibrium argument (from the Gibbs distribution). 

The length L present in the transition rate /x (Eq. 1281) or in the duration t is that of the potential plateau. 
It is the length of the domain minus the size of the two boundary layers on both ends L — 2SL. An analytical 
estimate of 6 L can be performed by calculating i/{dV/dy)y^o, where y is the position of the front, using the 
ansatz of equation |36] (see ? lA.ll) and letting L —>■ oo. The calculation is very similar to that of the potential 
difference AV. This gives a value of 8, while a more precise numerical calculation shows that SL ~ 5. 

3.3.4 Number of fronts 

At finite-temperature, the system may have several fronts separating the domain into different parts. In order 
to calculate how many, a useful analogy is that of a gas of weakly interacting subdomain fronts [35] ■ We 
can consider the probability Pn of having n fronts in the domain relatively to the probability po to be in one 
of the minima A^. We assume that the fronts do not interact at all. This assumption is realistic since in 
practice the interaction amplitude decreases exponentially with the distance between the fronts [T] |36] . The 
probability of having n fronts relative to po reads 

pL pL p L 

Pn = 2po dXi dX 2 --- dA„exp(-/3AK). (30) 

Jxi=0 Jx2=Xi J X„=X„-1 

Here the factor 2 accounts for the natural symmetry A —>■ —A and the integration is over all possible positions 
Ai, • • • , A„ of the fronts, the term AI4, = nAVi (see Ea fT2l) and is constant, i.e. independent of the front 
positions. The calculation of the above integral gives at the end 

jjn 

Pn = 2po—j-exp(-n/3AVi). (31) 

n\ 

Note that the factorial term 1/n! differs from earlier predictions m- It means that fronts behave in fact like 
identical particles in this treatment and it is often found in classical statistical physics. We now identify the 
regions of the plane (L, /3) for which observing n random-walk fronts is the more likely. It is simply done by 
identifying the curves = Pn-i (when having n fronts becomes more probable than n — 1) and = Pn+i 
(when having n + 1 fronts becomes more probable than n). We find 

Lexp(—/3AEi) — 1 < n < Lexp(—/3AVi). (32) 

It gives the size of the coherent domain L/n which is fixed by the competition between the level of noise 
present in the system and the potential cost of creating fronts. 

3.4 Theoretical phase diagram for reactive trajectories 

We provide here a compilation of the results we found previously. FigurelHlshows the type of reactive trajectory 
as a function of the inverse temperature /3 and domain size L. The different curves in this parameter space 
are found from conditions (I20L(I33|) as well as the results of FigjU Note that condition (1301) being exponential 
in /3, the existence of genuine instantons is concentrated in a very narrow band of relatively small size L, 
unless P is very large. Similarly, condition (1321) indicates that reactive trajectories with several fronts exist 
in a very narrow band of inverse temperature /3. Note also that in practice, most of the parameter plane is 
occupied by random walk type of trajectories. 


4 Numerical results 

We now present systematic computations of reactive trajectories using AMS algorithm. We discuss first the 
qualitative properties of these trajectories in S 14.11 We then discuss in a more quantitative way the trajectory 
properties and compare them with the asymptotic theoretical predictions. We decompose this discussion into 
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Figure 6: Sketch of the phase diagram of possible reactive trajectories in the {f3,L) plane. A random walk 
solution with one or two front is shown for illustration. 


three parts: mean first passage time (§ distribution of the duration of reactive trajectories (§ ig and 
escape probability IS 14.41) . The numerical procedure is described in the appendices IA.2I and IA.3I Along with 
the details on the algorithm, we explain the importance of the number of trajectories N and the number of 
realisations of the procedure on convergence in these sections. 

4.1 Reactive trajectories 

We consider here four different domain sizes L. We use N = 1000 trajectories in these simulations: this 
number is small enough for the simulations to be quick, but large enough for the outputs to be reliable. 

For L = 5 and /3 = 300, the transition rate 1/T is estimated from numerical results (See Appendix S IA.3.21 
Eg. [5T|) . In that case, it is l/T ~ 3.7 • 10“"^®. Note that in this system, the relaxation time is of order 1: this 
is much smaller than the mean first passage time that we computed. A typical reactive trajectory is shown 
in Fig. [7^. It corresponds to the first type of (flip) instanton going through the saddle Ag = 0 due to the 
relatively small size L. It is qualitatively comparable with the flip instanton shown in Fig. [5^). 

Another regime is for larger L. For L = 10 and /3 = 150 the transition rate is 1/T ~ 3.6 • 10“®^. In this 
case, one finds reactive trajectories with only one front (see Fig. [^o). It is consistent with the theoretical 
results where Ai becomes the saddle with lowest potential (see S 13.11 and Fig. Hb). The front propagates 
quickly from a; = 0 to a; = T/2 spends a very long time there and then moves forward to x = L. The main 
difference with a genuine front instanton (Fig. IJb) is that the position of the front fluctuates slightly around 
X = T/2. For these order of magnitude of the rate of escape, it is clear that a numerical approach such as 
AMS is necessary. Indeed, for spacially extended systems, direct numerical simulations are conceivable up to 
durations t < 10®. Computing mean first passage time that are much larger than that is not feasible. 

We then consider a domain of size T = 30 with /3 = 30 which yields 1/T ~ 8.3-10“^®. The algorithm finds 
reactive trajectories with one front as shown in Fig. [7](c). At this rather low /3, the solutions present some 
significative differences with front instantons (Fig. [3})). Although the transients from x = L to x ^ Li = 2b 
and that from X'^L 2 = btox = Q are very short as expected from deterministic relaxations, the front 
propagation has large fluctuations in the interval [T 2 ,Ti]. This is reminiscent of the three-stage process 
described in section 13.3.21 Note that this provides us an example of the potential boundary layer thickness 
(5T ~ 5 introduced in section [3.3.31 At this low /3 regime, saddle-point approximations of Freidlin-Wentzell 
and Eyring-Kramers theories are no more valid. In fact, most of the reactive trajectories found by the 
algorithm are not instantons when f3 becomes small. A typical case is shown in EiglTJi for /3 = 5 and T = 10 
where 1/T ~ 3 • 10“®. One still observes one front solution but with higly random motion. Another case at 
(3 = 7 and T = 30 (1/T ~ 4.3 • 10“"^), exhibits two fronts with rather regular propagation (see Fig[7^). This is 
comparable with trajectories minimizing the action at small duration T constraint [32j . Note that even at low 
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a) b) c) 



d) e) f) 


Figure 7: Examples of reactive trajectories computed using AMS algorithm for N = 1000, dt = 10 ^ and 
dx = 1/18. All the trajectories are typical of the average duration, a): L = 5, /3 = 300. b): L = 10, (3 = 150. 
c): L = 30, /3 = 30. d): L = 10, /3 = 5. e): /3 = 7, L = 30. f): j3 = 7,L = 110. 


/3, these two fronts trajectories are only a minor part of the fastest trajectories. Finally, we consider the case 
of a very large domain with L — 110 and (3 = 7 {1/T ~ 1.3 • 10“®). In this case, one observes several isolated 
fronts with random walk propagation until the system reaches Aq . Note that this situation only occurs for 
/3 < 7 and L >30. These values are critical in the sense that for f3 > 7 and L <30, reactive trajectories with 
a single front reappear. 

These first simulations show us that instantons are observed for large /3. The other expected diffusive 
regimes also found. Most of the relevant sets of parameters lead to transition rate so small that they cannot 
be attained using direct numerical simulations: methods such as adaptive multilevel splitting are of great 
help in such studies. 

4.2 Mean first passage time 

An efficient way to test the validity of our results is to compare quantities which can be precisely predicted. In 
particular, the mean first passage time gives us a good benchmark (see Eq. (|13p and the algorithmic formula 
used in Eq. (|5T|l . § IA.3.211 . We thus compute T as a function of (3 and L for /3 G [1,60] and L G [8,17]. In 
this regime, the saddle is Ai. We perform 10 independent realisations of the algorithm for N = 1000 and 
A^ = 10 000 in the case L = 17 and perform an average over these realisations. The results are shown in 
figure [8] (a). We superimpose the corresponding analytical predictions described in section [3.21 In practice, 
the numerical results given by AMS are estimates of the mean first passage time, trajectory duration, etc. As 
detailed in appendices IA.2I and IA.31 the actual values lie in an interval of confidence around this estimate. 
Using the mathematical convergence results mmaiiij and propagation of errors, we construct such intervals 
of confidence for the mean first passage times we computed and add them in the form of error bars for figure [5] 
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Figure 8: (a): Logarithm of the mean first passage time computed using the asymptotic formula (EqlTS]) and 
using the AMS one fEa lSTj) as a function of /3. (b): Logarithm of the mean first passage time as a function of 
/3 zoomed in the interval 9 < /3 < 26. (c): Logarithm of the mean first passage time (analytical and numerical 
result) minus AF/3 as a function of ln(L — 25L) for several values of 20 < /3 < 30. 


(a). We took care to use enough trajectories N and enough realisations of the algorithm for these interval of 
confidence to be small, explaining why error bars are not visible in the figure. In the rest of the article, these 
intervals of confidence will not be included for readability. 

We move to the discussion of the physical content of the numerical results. We can see in figure |S] (a) and 
more particularly in the zoom of figure |8](b) that for /3 > 7 and L < 13, we have a near exact match between 
the mean hrst passage time estimated by AMS and the analytical result of the Eyring-Kramers formula. The 
ratio of the numerical estimate over the analytical prediction is 1 ± 0.1. This means that in this part of the 
parameter plane, the transition rate 1/T decays exponentially with j3 & rate — AF and also exponentially 
with L at a rate —1/-\/2. The Eyring-Kramers formula is an asymptotic result in the large /3 limit, so that 
the small discrepancy at small (3 is not surprising (this will be more visible when we examine the results on 
the escape probability in section 14.41) . In section 13.31 we explained why the Eyring-Kramers is not valid for 
larger sizes. We expected that the agreement between numerical estimations of the mean first passage time 
and the prediction of the Eyring-Kramers formula would only be valid up to a certain size. The value L = 13 
gives us a lower bound on the length L*(/3) above which these large /3 result should not be expected. This 
lower bound is in practice independent on /3: this is in perfect agreement with the prediction that L* grows 
logarithmically with (3 (see Eq. (l20l) L 

For L > 15, the numerical estimate of the mean first passage time is smaller than the value given by 
the Eyring-Kramers formula as shown in figure |8] (b). This indicates that L has crossed L* for this range 
of p. Nevertheless, we find that the Freidlin-Wentzell large-deviation prediction is very robust although the 
Eyring-Kramer prefactor is no more valid m- The logarithm of the mean first passage time still grows 
linearly with /?, with the same slope AV as the asymptotic Eyring-Kramers result. Meanwhile the ratio of 
the numerical estimate over the analytical result is approximately 0.5 for L = 15 and approximately 0.4 for 
L = 17. This value is independent on /3. This numerical result is reminiscent of the crossover between an 
exponential growth of T with L (see Eq. (1281) § 13.2.21) and a much slower polynomial growth of T with L (see 
Eq. (UHl), ^ 13.3.21) 

We examine the large size L regime in more detail in view of the analytical results of section [3.3.2l fEa. dT51) L 
In the diffusive regime, this formula tells us that ln(T) — PAV = ln(L — 26L) + ln(/3Zi/16). The boundary 
length 6 L had been introduced in section [3.3.31 to account for the actual length of the potential plateau. We 
display E = ln(T) — PAV as a function of ln(L — 2SL) in figure [8] (c), where T is estimated using the AMS 
and AV is computed numerically. Three values of ,5 = 20, 25 and 30 were considered. We can first see that 
E is independent on /3, confirming that T grows exponentially with /3, whatever the size. We find an affine 
growth of E with ln(L — 26L). This indicates that in the large size L regime, the mean first passage time 
T grows as a power of L. This finding differs quite strongly from the exponential growth expected from the 
Eyring-Kramers formula. A fit of the numerical data gives an estimate of 1.5 for the slope of A as a function 
of ln{L-2SL). 
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Figure 9: (a): Distributions of normalised durations of reactive trajectories for L = 10 in a range of /3 
compared to a Gumbel distribution, (b): Average duration and variance of reactive trajectories as a function 
of ln(/3) for L = 8. (c): Average duration of reactive trajectories as a function of ln(/3) for 8 < L < 15. (d) : 
Variance of duration of reactive trajectories as a function of ln(/3) for 8 < L < 15. 

4.3 Duration of reactive trajectories 

We now investigate the statistical properties of the duration of reactive trajectories. As explained in sec¬ 
tion 13^531 the duration of a reactive trajectory is the time elapsed between the instant where the trajectory 
leaves C and the instant where it reaches dB (Fig. [1]). A realisation of the AMS with N trajectories will thus 
give us N durations. From theses N durations, we can construct an empirical distribution of the trajectory 
durations (Fig. [31), estimate the mean duration r (Fig. [3] (b,c), Fig. (TU] (a)) as well as the variance of the 
distribution (Fig. [3] (b,d). Fig. [13] (b)). In practice, the quality of the sampling of the empirical PDF as well 
as the estimation of the mean and variance depends on N. One can verify numerically that we have central 
limit type interval of confidence on the cumulants m- We chose N = 10000 in the systematic study so that 
this interval of confidence is small enough. 

We first focus on the non diffusive regime, at large (3 and small L. We compare the numerical results 
to the analytical results for transition over a saddle in a one of freedom system discussed in section 13.2.31 
Figure [3] (a) shows that the distribution of normalized durations is independent of /?. The comparison with 
a normalized Gumbel distribution (ITS)) is relatively good good, in particular for t > t. We notice a slight 
systematic difference for t < t. In order to go further with the comparison with the one degree of freedom 
case, we displays the average duration of trajectories and the variance of the distribution as a function of 
ln(/3) in figure [3] (b) for L = 8 . The data are consistent with r being an affine function of ln(/3) and the 
variance of the distribution being nearly a constant, as in the one degree of freedom case. The estimate of the 
slope of t(/3) is 4.7 while is variance is 5.5 ± 0.2. These values fall within 10% of the one degree of freedom 
prediction of a slope of t(/3) being l/|As| ~ 5.3 and a variance being tt/(|| v^) ~ 6.8. 

However, as the size L is slightly increased, the average duration and the variance of the distribution 
departs much quicker from the asymptotic /3 —>■ oo regime than the mean first passage time. The average 
duration as a function of ln(/3) is displayed in figure [3] (c) for 8 < L < 15. One notice that for L = 15, we 
are already leaving the regime of r affine in ln(,5). Moreover, the slope of r Vs. ln(/3) becomes clearly larger 
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Figure 10: (a): Average duration of the reactive trajectories rescaled by {L — 25L)"^ as a function of /3 for 
domain sizes larger than 30. (b) Variance of the distribution of reactive trajectories durations rescaled by 
and (L - 25Lf as a function of {3 for domain sizes larger than 30. The additional blue dots are an 
example of rescaled average duration. 


than 1/1 As I for L = 10 and L = 13. The departure is even more striking when we consider the variance of the 
distributions of durations as a function of ln(/3) for 8 < L < 15 (Fig. [9] (d)). While we do find the constant 
variance regime for L = 8, we notice that for L > 10, we find a variance which is nearly affine in ln(/3), so 
that var oc r. This situation is intermediate between the diffusive and the non diffusive regime (see S 13.3.31) . 
This illustrates the non-triviality of the question of the distribution of durations of reactive trajectories in 
stochastic partial differential equations. 

We then move to the fully diffusive regime. It occurs for domain sizes larger than L > 30 in the range of 
P available to numerical simulations. Like in the non diffusive regime, we compute the average and variance 
of the duration of reactive trajectories (Fig. (TUI). As explained in section 15.3.31 we expect trajectories to be 
random walks on the potential plateau and the trajectory durations to have the corresponding behaviour. 
Computing r rescaled by {L — 25L)'^ as a function of /3 for 30 < L < 100 strikingly confirm the result (Fig. [TU] 
(a)): the values for L > 45 all fall on a master line. Note however that the slope is nearly 1/6, as would be 
expected if the diffusion coefficient were 1//3 and not 1/16, as would be expected if the diffusion coefficient 
were 3/(8/3) [37]. The discrepancy with the mathematical result remains to be explained. The variance of 
the trajectory durations is a growing function of /3 which again corresponds very well to the random walk 
behaviour of the front (Fig. (TOj (b)). Indeed, when rescaling it with (L — 25L)^, it is linear in /3. Moreover a 
second rescaling by \f5j^ shows remarkable agreement with the average duration, indicating that var/r has 
the expected ratio. 

4.4 Escape probability 

We now consider estimates of the escape probability a (see also Appendix IA.2I) . This is the probability to 
reach the set B before the set A starting from the iso surface C (see Fig. [H (a)). It can formally be written 
with P (see Eq. lU § 13.1|) by integrating over all the initial conditions uq on C and all the final conditions 
iti on dB [21j . Unlike the transition rate computed in the previous subsection, this quantity depends on the 
definition of the sets A, B and C. The escape probability has less physical importance than 1/T. However, 
a corresponds to the committor (see [201122] 1 whose asymptotic computation is difficult in many degrees of 
freedom systems. Since the committor is the solution of the same type of backward Fokker-Plank equation 
as the mean first passage time, large deviations results exist showing that it should decay exponentially like 
exp(—/3AU) [5TJ[ni], with AU the potential difference between the starting point on C and the lowest saddle. 
Moreover, this quantity is of fundamental numerical importance, as it is required to compute the transition 
rate. Besides, it is the only quantity on which mathematical results on the convergence of the algorithm exist 
(see [m |44| ESI |46| |47l [48] and the appendix). 

We study how a depends on /3 for sizes L G [4.5,30]. The number of trajectories is A = 20 000 and 
the number of independent realisations of the algorithm is taken equal to 20. Figure [TT] (a) clearly shows a 
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Figure 11: Examples of Ina as a function of /3 with the corresponding linear fit: (a): L = 7.07. (b): L = 17.3. 
(c): as a function of L. (d): Ordinate at the origin of ln(Q!) at /3 = 0 as a function of L. 


change of regime with at /3 ~ 8 . Above this value, we observe a large-deviation result where In a decreases 
linearly with /3. We investigate more in details this behavior by dividing the corresponding slope by AV = 
V (Asaddie) - F(C) ~ F(Asaddie) “ V (A^) using either A^addie = Aq or Agaddie = Ai depending on the values of 
L. This is shown in figure ITT] f cl. The quantities are very close to the value -|-1 at ±5% precision which agrees 
very well with the large-deviation predictions in these large /3 regimes. Note that these results are obtained 
for ,5 > /3* nearly independently of L. A large deviation principle exists for the escape probability (Fig. [TT] 
(a,b)), while reactive trajectories are no more instantons (see FigslTt and[3}D). Once again, it means that the 
escape probability is essentially related to the cost for creating one front and it is a very robust mechanism 

132]- 

A significant advantage of AMS algorithm is that it gives precise information on Eyring-Kramers pref¬ 
actors. They can be obtained by computing the ordinate at the origin of In a as a function of the size L 
for instance. This is shown in Fig. [TT] (d) with prefactors belonging to [0.1,10]. The variation with size of 
the ordinate at origin is likely related to the numerical precision on the slope: an overestimate of the slope 
in absolute value leads to an overestimate of the ordinate and vice versa. However at hxed /3, one generally 
notes a decrease of a with L. 

4.5 Phase diagram 

We give a synthesis of these results by constructing a numerical version of the diagram shown in Fig. [51 One 
can determine the number of fronts through visualisation of the space-time representations. The discrimina¬ 
tion between a flip and a one-front trajectory is more easily done by looking at the average duration (see Fig. 
ist). However, visualisation alone is not enough to reasonably determine whether a trajectory is an instanton 
or a random-walk trajectory. For a given point (/3, L) in parameter space, we consider a trajectory to be 
an instanton if the average duration and the mean first passage time falls into the non-diffusive regime (see 
Fig. [5}: Vs. Fig. [TUI and Fig. |5|). For most sizes, the transitions belong clearly to one regime or another. We 
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Figure 12: Phase diagram of possible reactive trajectories in the (/3, L) plane using AMS algorithm to be 
compared with figure [S] In particular, we use the same notations: Iq: instanton flip, Ii: instanton with one 
front, Rq: random walk flip, Ri: random walk with one front, R 2 : random walk with two fronts. 


assumed that the transition was out of the low noise regime if the relative difference between the mean first 
passage time calculated using AMS and using the Eyring-Kramers formula was more that 15%. The precise 
position of the boundary between regime will depend on the value of this tolerance. This criterion gives a 
lower bound ,5 > 12 and an upper bound L < 13 for the instanton regime. All these informations are encoded 
in the phase diagram of figure fT^ There is a very good agreement with the theoretical predictions both on the 
different regions of parameter space and the type of reactive trajectories (Fig. IS]). In our simulations, we were 
not able to obtain reactive trajectories with more than two fronts. The exponential curves at the boundary 
of some of the regions are not resolved here and appears to behave as straight lines, as a consequence of the 
rather poor grid resolution of the (/3, L) plane. 


5 Conclusion 

This article intends to present a complete theoretical and numerical analysis of transition rates and reactive 
trajectories between two metastable states of the 1-D Allen-Cahn equation. In this work, we have used a 
recent algorithm called adaptive multilevel splitting dedicated to the computation of very-low probability 
events. We are able to numerically investigate metastability of the Allen-Cahn equation in a fast and efficient 
way by systematically computing reactive trajectories using this algorithm. This is one of the main aspect 
of this work. It allows us to synthesise existing asymptotic theories like Freidlin-Wentzell large deviations 
m, Eyring-Kramers law of mean first passage time [Ml HH121 HD, and theory of front motion [21137]. The 
numerical results validate these theories but also bridge the gap between them and go further. 

The combined use of numerics and theory shows that in the zero noise limit, reactive trajectories are either 
global reversal of the flow when the domain size L is small, or a front nucleation at one end of the domain 
and its propagation toward the other end provided L is sufficiently large. Solutions with several fronts are not 
global minima of the action |32M36j . and are never observed numerically for large enough /3. We are also able 
to compute with good accuracy both the escape probability and the transition rate in many different regimes. 
We show that the escape probability has a very robust large-deviation behavior in the limit of zero noise, 
as expected from Ereidlin-Wentzell theory. Moreover, in regimes where Eyring-Kramers law does not apply 
anymore, we always observe a large-deviation behavior of the transition rate whenever the escape probability 
does. It calls for possible further investigations on the prefactor alone in these situations. Another remarkable 
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result is the robustness of the large-deviation behavior of the escape probability at finite-amplitude noise and 
when the most probable trajectories are no more instantons. Theoretical predictions, especially those on size 
dependence, are in fact valid in a much larger range of parameters than the low noise limit /3 —>■ c». The 
validity of the prediction merely requires the creation of a single front and a long enough reactive trajectory 

m- 

The AMS algorithm is also able to provide accurate statistics on the average duration r of reactive 
trajectories for which no theory exists yet. Remarkably, we find similar scaling laws in Allen-Cahn equation 
than for systems with only one degree of freedom, for which theoretical analysis exist [50]. We obtain that 
r oc ln/3 X exp(L/-\/2) when the curvature at the saddle is large, and r oc /3L^ when it is small. 

This work supports the idea that AMS is a very efficient tool for analysing extremely rare transitions in 
systems with a large number of degree of freedom and SPDEs as well. It is a promising approach for studying 
more realistic systems arising from practical problems lii. Although Allen-Cahn equation has a rather 
simple structure as being a gradient system, AMS is likely to give useful insights for non-gradient systems 
as well, as there are no principle limitations in its use. In particular, an interesting step is to consider 2-D 
non-gradient SPDEs such as the stochastic 2-D Navier-Stokes equations |4] whose theoretical description is 
much more subtle than for gradient systems [30]. We are confident that AMS approach will becomes more 
popular in the future due to its simplicity, robustness and efficiency. 


A Appendices 

We first develop the calculation of the approximation of the negative eigenvalue of the Hessian of the potential 
V We then give here some details on the AMS algorithm, first we provide a detailed description of the algorithm 
itself, then some basic applications like the computation of mean first passage time. We end with a practical 
discussion on some numerical convergence issues. 


A.l Approximating the Eigenvalue 

In this first appendix, we detail the calculation of an approximation of the first eigenvalue As of the Hessian of 
V, presented in section [3.2.21 This eigenvalue controls the prefactor of the mean hrst passage time (Eq. [14]). 
The eigenvalues Xi and eigenmodes are the solution of the problem 


-52$, + ( 3 A 2 0 - 1) $^ = X^,A,,o^^,As 




(33) 


with boundary conditions $^(0) = ^i{L) = 0. The eigenvalue problem can be posed either at the saddle, for 
As, or in the minimum Ag. The operator i7s,o = 0 hessian of the potential V. It is self adjoint 

on the set of function (regular enough) that cancel out on both ends of the domain x = 0,L. It is well known 
that in the limit of infinite size L —>■ 00 (see |36j and references within): we have 


As = Ai_oo = tanh 



$00 = 1 ~ tanh^ 




Aoo =0. 


(34) 


This eigenmode is a goldstone mode, related to the translational invariance of the front in an infinite domain 
[55] . We denote hyperbolic tangent 'Hu[x) = tanh((a; — u)/\/2), which will simplify our notations. 

How As converges toward 0 in a finite domain of size L ? In order to address that question, we will provide 
an upper bound for Ag. For that matter, we use the fact that Hs is self-adjoint. Its Rayleigh quotient 

R{Hs,f) = , /(O) = fiL) = 0 , (35) 

Jo f 

then verifies the min-max theorem, which states that R{Hs,f) is bounded by the lowest (As) and highest 
eigenvalue of Hs. Moreover, we have that R{Hs, $) = As and that given an approximation of $, R{Hs, 4') 
will converge (from above) toward As as dt converges toward 4>. 
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We will now tailor an analytical approximation of $ and Ai so as to give an upper bound on As- This 
upper bound will serve as an approximation. The lowest eigenvalue As and the eigenmode $ correspond to 
the unstable direction of V at the saddle: small motion of the position of the front near L/2. If the field A 
is slightly displaced away from the saddle Ai along the unstable direction, the front (located at position y) is 
displaced by a small 6 x = y — L/2. We then have A^ Ai + (5a:$. In order to obtain an approximation of 
we start from an analytical approximation K{x,y) of the field ^(a;) when a front is located at position y 


K{x,y) = 'Ho{x)'HL{x)'Hy{x). (36) 

Performing the expansion gives us 

« ^ (a:, Y + fe) = K + 0{Sx'^). (37) 

This leads us to choose 

^ = -V2^ ^a:, = -'Ho{x)'Hl{x) (^1 - . (38) 

The minus sign and the square root of 2 simply ensure that 'h is a positive function whose maximum is I at 
X = Ll2. In order to give a good analytical approximation of the ratio of integrals I = R{Hs, d>), we set 

- So “ 1 )) 

A, = no{x) + nL{x)-nL{x),i= -^— -j-— -^-. ( 39 ) 

/o 


Note that $oo does not cancel out at x = 0,L, nor is Ai strictly the saddle point of V. As a consequence 
/ is not strictly speaking a Rayleigh quotient, nor should we a priori expect to verify exactly its properties. 
The ratio / is a good approximation of I because like decrease like 1/cosh'^ away from T/2, up to 
corrections on the boundaries (Fig. 1131) . Any relative contribution to the integrals near x = 0 and x = L, 
where the errors are made, are exponentially small by a factor exp"‘(—L/(2-\/2)) = exp(—Lv^) relatively 
to the contributions from the neighbourhood of x = L/2. Meanwhile, the relative difference between $oo 
and 'k, or between Ai and Ai are also exponentially small near L/2. Note however, that it is fundamental 
that Ai has the proper boundary conditions at x = 0 and x = L, as it is they that will contribute to the 
exponential decay of the eigenvalue. The quality of this approximation depends on how large L is. For L = 10, 
exp(—L-\/2) — 7 • the system is quickly in the range of validity. 

The denominator is calculated in the same way as AVi, it is equal to (4-\/2)/3. In order to calculate the 
numerator, denoted J, we first note that 


i7$oo =- 


We process the sum 


3(i-Hl(x)2) \no{x) +'Hl{x)\ I -'('Ho(x) +-HL(x))+2'Hr (x) | . 


=s 


=s 


S = 2e ^ 




1 + exp (•\/2(x — L)) + exp (—v^x) + exp {—^/2L) 


(40) 


(41) 


In the neighbourhood of T/2, only 1 has a significant contribution in the denominator, so that this sum can 
be approximated by the numerator in H^oc- Note that it is negligible compared to tanh((x — L/2)/-\/2): 
Both cancel out in L/2, while only the hyperbolic tangent is of order 1 in the direct neighbourhood. Now J 
reads 


J = 


-12e ^ 




'Hl (x) 

2 ^ ^ 



2x-L\ 

J 


— exp 


2x- L\\ 

JJ 


dx 


(42) 
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Figure 13: Logarithm of the eigenmode and its approximations in a domain of size L = 20. 

We eventually note that the only contributions to the integral are around L/2, so that we can let the bounds 
go to infinity. This gives us J = —32-\/2e“^/'^, then 

/ = -24e"^ . (43) 

This is an approximation at first order in exp(—L/-\/2): the corrections are exponentially small as L takes 
large values that can however be considered in AMS simulations. A numerical calculation of the integral in 
/ confirms this assertion. 

A.2 Algorithm description 

The general setting is to have a Markov process describing our system, i.e. (Af“)t>o with = xq & £■ The 
phase space 8 of the system can be either finite-dimensional of infinite-dimensional. We assume that there 
are two sets A and B with AC\B = %. The goal is to estimate the probability a to reach B starting from the 
initial condition xq before going to A. If tc denote the stopping time defined as rc(x) = inf{f > 0, Xf S C,} 
for a given Borel set C, the problem translates into computing 

a = a(xo) = Pr(r£5(xo) < r^(xo)). (44) 

A reactive trajectory is a particular realisation of (1441) . Note that this probability, when seen as a function of 
xq, takes values zero ii xq € A and one if xq G B. It is called committor or importance function or equilibrium 
potential in mathematics. In particular, for diffusive processes, it solves the backward Fokker-Planck equation 
with boundary conditions a{x) = 0,x S A and a(x) = l,x € B. The algorithm not only gives a pointwise 
estimate of the backward Fokker-Planck but also provides the ensemble associated with the probability. It 
is clear that for high-dimensional systems, access to the Fokker-Planck equation or doing direct Monte-Carlo 
simulations when a is too small is, most of the time, out of question. 

The main idea is to introduce a scalar quantity which measures how far a trajectory is escaping from the 
set A before returning to A. There is no unique way to measure how far a trajectory is escaping from the 
return set A, or put in another way there is no unique order relation in the space of trajectories. However, a 
natural choice is to consider the following quantity 

Q(x) = sup $(Af), $ : E ^ R. (45) 

The function $ which is often renormalized to be in [0,1] is called reaction coordinate or observable and is 
such that 4>(A) = 0 and ^{B) = 1. Note that the choice of $ might become crucial and we will comment it 
hereafter. The algorithm is doing the following (see also Fig. [T]): 

Pseudo-code. Let iV be a fixed number, and $ given. 

Initialization . 

Set the counter /C = 0. Draw N i.i.d. trajectories indexed by 1 < f < A^, all starting at xq until they 
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either reach A or B. Compute Qi(xo) for 1 < z < iV. In the following we will write Qi these quantities 
since their initial restarting conditions might change. 

DO WHILE fminQ, < b) 


Selection', find 


Let Q = Q^*. 


i* = argmim Qi. 


Branching', select an index / uniformly in 
{I, • • • TV} \ {z*}. Find t* such that 


C = inf{$({Xaz) >Q} 

Compute the new trajectory z* with initial condition x* = {Xt*}i until it reaches either A or B. 
Compute the new value Qi*{x*). 

/C^/C + 1 


END WHILE 

The performance of the algorithm depends rather crucially on the choice of $. In fact, it is possible to show 
that when $(a:) = a{x) (see eq.(|44])) or when the system has only one degree of freedom {£ = R), the number 
of iterations /C when the algorithm has stopped has a Poissonian distribution [46l [45] 

/C ~ Poisson(—TV In a). (46) 

In other words, the quantity yields an estimate of — In a (see m for particular SDEs). Moreover, 

(47) 

gives an estimate of a itself. In some cases, one can demonstrate a central limit theorem for a, with a variance 
scaling like 1 / 

\ Oi J N—¥oo 

Some sets of hypotheses lead to a bias behaving as smiE]. With more loose hypotheses or in more complex 
cases, the variance still scales like 1/'/N, but can be larger than what is predicted m- 

Note that for arbitrary $ and in dimension larger than 1, the statistical behavior of the algorithm is a 
difficult mathematical question and it is still open [45] . In particular, it is unclear whether these estimates 
are biased or not, even in the limit of TV —>■ oo. In practice, a “bad” reaction coordinate which differs 
significantly from the optimal committor, has a fat-tailed distribution although its variance still behaves like 
Therefore, the larger TV, the better precision one obtains. Overall, the algorithm numerically exhibits 
very robust results. There are other formulation of the algorithm, in particular, it is possible to kill a given 
proportion of the TV trajectories at each algorithmic step (typically, if one kills p trajectories at each step 
followed by replacement, an estimate of the probability becomes (1 — ;^)^). This version has the advantage 
to be computationally more efficient when parallelized f|47jL 

A.3 Practical applications for Allen—Cahn equation 
A.3.1 Reaction coordinates 

In our case the set A is a neighborhood of the stationary solution A(} and the set B a neighborhood of Aq . 
The choice of the reaction coordinate is often dictated by the problem itself, a natural one is to choose a 


25 






quantity which tells how far one is from . We use the following reaction coordinates 





$„(A) = ^i{A) 


1 < A,Ai> 

7 l Pill 


( 49 ) 


where A = A(x) dx and < A, B >= A{x)B{x)dx,\\A\\'^ = A^{x)dx. These reaction coordinates 

have the property that they vanish at Aq and are equal to 1 at A^ . The second reaction coordinate gives 
some additional weights to the fronts present in A. In practice, we chose 


dA = = 0} , = {A\$„(A) = 1} , C = {A\^»(4l) = 0.05} . (50) 

For domain of large size {L — 60), both values = 0.05 and = 0.025 were tested. The two definitions led 
to differences in computed quantities (3, r) clearly smaller than their variance. For this reason, we consider 
than this change of definition of C does not play any role on the final result. However, this change plays a 
clear role in the computation time of the generation of initial conditions distributed on C, especially when /3 
or L are large. For this reason, = 0.025 was systematically used for computations at L = 80 and L = 100. 


A.3.2 Statistics from the output 

There are two important quantities which can be obtained from AMS output: 

• distribution of duration of reactive trajectories. In ST], a numerical analysis shows the convergence of 
this distribution. In the case of SPDE and Allen-Cahn for example, there is no theoretical results on 
this distribution. In our case, one can check for instance the rate of growth of these lengths with /3. 
More generally, since AMS is providing an ensemble of reactive trajectories, it is possible to perform 
many a-posteriori statistics which become accurate as N is larger. 

• mean first passage time In order to compute this quantity, one defines a isosurface C surrounding the 
set A. One generates an ensemble of N trajectories which all start at A^ . One then estimates the mean 
time Ti it takes to cross the surface C and T 2 the mean duration of nonreactive trajectories (that is, the 
length of trajectories between the hitting time of C and the return time to A) [62] . Assuming that there 
are a mean number < n > of nonreactive trajectories, the mean first passage time T is 

T =< n > (ri + T 2 ) + (ti + t^), (51) 

where the mean duration of reactive trajectories found by the algorithm. The number of failed 
attempts < n > is in fact related to a by a = . 

A.3.3 Convergence issues 

There are two sets of parameters to take into account: first, the SPDE discretisation with timestep dt and 
spatial resolution dx, second, the number N of the trajectories in the algorithm and the choice of the reaction 
coordinate $. We first illustrate the convergence properties of the SPDE itself. We use a test case of L = 7.07 
using N = 20 000 trajectories and the norm reaction coordinate (which yields better results than [47]). 
We increase the spatial resolution from dx = ^ to dx = -^ and decrease the timestep dt accordingly. Note 
that, in order to ensure the stability of the numerical integration of the diffusive part, the constraint is to 
have D > . The choice of dt in addition to this constraint crucially depends on the accuracy of hitting 

times (of A and isosurfaces of the reaction coordinate $) needed in the algorithm. It is well-known that these 
hitting times converge like \fdt |63] bringing a limit to the accuracy of the solution and the probability a 
found by the algorithm as well (see [47]). 

Estimates of In a iFiglTihl and the average duration of reactive trajectories (FigfHbl are obtained 
for the different parameter cases. The result is that increasing the spatial resolution does not improve the 
accuracy of In a whereas it improves the accuracy on the average duration of reactive trajectories (Fig. 
(Hb). Nevertheless, the qualitative behavior of is robust to the parameter dx and dt. In this article, we 
have decided to keep the values dx = ^ and dt = 10“^ in all the simulations. 
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Figure 14: (a): Logarithm of the escape probability as a function of /3 for L = 7.07, TV = 20 000, and 
three different values of dt and dx. (b): Average duration of reactive trajectories, in the same setting, as a 
function of ln(/3) for the same three resolutions, (c): Comparison of the rescaled variance uo of the escape 
probability as a function of the logarithm of the escape probability, for the two reaction coordinates and 
Both Allen-Cahn equation and the two-saddle model (see Eg. (1551) 1 are shown. 


We now discuss the effects of the number N of initial trajectories and the choice of the reaction coordinate 
$ in the quality of the results. In practice, for most choices of $, the results converge in the limit N —>■ oo. 
However, the bias and variance of the escape probability estimate can behave rather differently depending on 
/3 and the choice of $ [47]. We look here at the relative variance 


ao = cr 


Vn 

ay/— In a 


(52) 


as a function of /3 and for the two reaction coordinates and Here a is the variance of the estimate of a 
obtained by the algorithm. Note that in dimension 1 or when $ is the committor, we have ctq = 1 (mills]). 
In order to simplify a bit the discussion, we also perform an analysis on the following system with two degrees 
of freedom. 


d{x, y) = -VVdt + ^{dWudW 2 ),Vix, y) = ^ ^ - 0.3 ^ + x^y^^ . (53) 

The main interest of this reduced model is that it displays a similar behavior of ao as the escape probability 
a becomes small. The reaction coordinates for (|53|) are ^i{x,y) = and $„(a;, 2 /) = \yj{x + lY -\- y'^ 12. 
We used the same definitions of A, B and C as equation. (ISOll . Note that this constitutes a slight change in 
the dehnition of A and B compared to what was used in [47]. Figure [HJ: shows the behavior of In a as a 
function of (3 for both Allen-Cahn equation and the reduced model. 

For the linear reaction coordinate the variance of the estimate distribution increases with /3 and has 
a bias going to zero when N ^ oo. This bias increases when /3 —^ oo. For the Euclidian reaction coordinate 
one obtains better results in the small model with a smaller variance and bias. These results appear to 
be generic among models m and is the subject of further theoretical investigations. Based on these results, 
we perform all the Allen-Cahn simulations using the reaction coordinate 4>„ in (|49p . These are done with a 
rather small number of trajectories N = 1000, when displaying qualitative properties, and N = 20 000 for 
more quantitative results. 
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